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A Monte Carlo scheme to sample the screening potential H(r) of Yukawa plasmas notably at 
C*| , short distances is presented. This scheme is based on an importance sampling technique. Com- 

parisons with former results for the Coulombic one-component plasma are given. Our Monte Carlo 
simulations yield an accurate estimate of H(r) as well for short range and long range interparticle 
£h , distances. 
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I. INTRODUCTION 



'-^J , In this work we present a Monte Carlo scheme devised to compute the pair distribution function g(r) of a strongly 
£h ' coupled plasma, notably at short distances, i.e. for values of r smaller or of the order of the ionic radius, with a high 
accuracy. The model considered in our study is the Yukawa One-Component plasma (YOCP), i.e. a system made of 
i N identical classical ions of charge Ze immersed in an uniform neutralizing background of electrons. The effective 

interaction between two ions is supposed to be of the form v a (r) — (Ze) 2 y a (r) where y a ( r ) = exp(— ar)/r (a > 0) 
t-H , is the Yukawa potential and a the screening parameter. In the limit a — ► the YOCP reduces to the well-known 
coulombic one-component plasma (OCP) 1 . The configurational potential energy in a domain A of the ordinary space 
Q\ ! M 3 can thus be written as 2 
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where pg = —NZe/A is the uniform charge density of the background. The constant £ which appears in the r.h.s of 
I ■ eq. (1.1) fixes the zero of energy and reads as 

§ ; £ = \h I^ r ) - ;1 = - a l 2 ■ ( L2 ) 

\ In the thermodynamic limit the thermodynamic, structural, and dynamical properties of the YOCP depend solely upon 
■ two dimensionless parameters, namely the coupling parameter T = f3(Ze) 2 /a and the reduced screening parameter 
a* = aa where (3 = l/fc^T (ks Boltzmann constant, T temperature) and a is the ionic radius (47r/oa 3 /3 = 1, where 
^ p = N/ A is the number density of particles). 

The thermodynamic properties of the YOCP are well-known nowadays thanks to extensive Monte Carlo (MC) 
simulations performed cither within periodical 3 or hyperspherical 2 boundary conditions. Much less is known about 
g(r) and the related screening function H (r) defined as 
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g(r) = exp (-Tay a (r)+TH(r)) 



(1.3) 



-ff(r) plays an important role in estimating the enhancement factors for the thermonuclear reaction rates 4 . As for 
r — > g(r) ~ exp(— ray Q (r)) the values of g(r) for r — > are extremely small for large T's, which precludes a numerical 
study by means of standard MC simulations and biased MC schemes are therefore unavoidable. Such a scheme is 
presented in next sec (II) , it is a synthesis of two biased schemes applied formerly in the determination of the cavity 
function of hard spheres 5 in the one hand, and in the calculation of the screening function H(r) of the OCP in the 
other hand 6 . 

Results of MC simulations are reported in sec (III) both for the OCP and the YOCP cases. They must be considered 
as preliminary results, our ultimate goal being to establish a complete data basis of screening functions H(r) for a 
wide range of (r, a*). 



II. SAMPLING THE SCREENING FUNCTION 

In the canonical ensemble the pair distribution function g(r) is given by 



J Il ti dn S(r- f 12 ) exp (-0V A (1, 2 ■ ■ ■ , N)) , 
JYlZidn exp(-/?14(l,2...,iV)) 



g(r) = A( J lii=1 r h „ 2 , * „J ) ■ (2-1) 



where the brackets < . . . > denote a canonical thermal average and A is the volume. In cq. (2.1) r 12 = r± — r 2 and 
we have implicitly assumed that the system was homogeneous which is verified if periodical boundary conditions are 
adopted. In practice g(r) can be computed in a standard MC calculation with a good precision only for r > r m i n . 
For instance at T ~ 100 we have typically r min /a ~ 1. In order to compute g(r) for r < r min we follow the suggestion 
of Ogata 6 and rewrite eq. (2.1) as 



= A J 11=1 dfj S(r- r 12 ) exp(-/?y A (l, 2 - ■ ■ ,N) - f3w(r 12 )) exp(0w(r 12 )) 

I dn exp; -i\v\. 2 • • • , AT)) 
« g w (ri 2 )exp(+(3w(r 12 )) . (2.2) 



In (2.2) w(r) is a priori an arbitrary function and Ogata has shown how to take advantage of that to devise an efficient 
MC biased scheme. 

The function g w supports the following simple physical interpretation. Let us consider a mixture made of (N — 2) 
Yukawa charges and two test particles labelled (1,2). These two particles interact with the (N — 2) other ions via 
Yukawa potentials but their mutual potential energy is defined as w(r^l2)) + v a (ri 2 ). For practical purposes, it is 
clearly clever to choose (3w(r) = —Tay a (r) + TH(r) where H(r) is a good estimate of the true screening function 
H(r). Indeed it follows from eqs. (1.3) and (2.2) that g w (r) oc exp (r [H (r) — H (r)] ) . As a consequence, if H ~ H, 
then g w (r) is practically constant and therefore easy to determine numerically. However we are only half the way since 
g w {r) is known only up to a multiplicative constant. This normalization constant can be reexpressed, as discussed 
by Ogata 6 , as a difference of free energies which can be determined in the course of a MC run as thermal averages. 
However, in the MC simulations of this author, the variance on these averages turn out to be quite large which 
precludes an accurate estimate. Therefore we adopted the method that Patey and Torrie devised for computing 
the cavity function of hard spheres 5 . Suppose that g w (r) has been computed by a MC simulation of the mixture 
described above for, let say, < r/a < 2. In practice it can be achieved by the choice (3w(r) = —Tay a (r) + TH(r) 
for < r/a < 2 and (3w{r) = 00 for r/a > 2. Then an unnormalized pair distribution g u (r) = g w {r) x exp(/3w(r)) 
can be computed in the range < r/a < 2 and compared to the normalized go(r) obtained by a standard MC 
simulation. This comparison in the range of distances where both g (r) and g u {r) can be determined precisely (i.e. 
Train/ a < r/a < 2) yields an accurate determination of the wanted normalization constant. Therefore both g(r) and 
H(r) can be computed accuratly for all r's. 



III. MONTE CARLO SIMULATIONS 



The scenario of sec. II was scrupulously applied in our MC simulations which were performed within hyperspherical 
boundary conditions 2 . For each pair of (T,a*) a standard and a biased MC simulation was performed. The runs 
were divided in typically ~ 20 sub-runs in order to compute the statistical errors by a method of blocks 7 . In the 
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preliminary data reported here the considered systems involved N = 500 particles but a systematic study of finite 
size effects is under way. 

In the case of the OCP we chose for H(r) the fit of H{r) provided by Ogata 6 . However a small attractive potential 
8H(r) was added to H(r) at short distances to enhance the sampling of g w (r) at small r. The form of 8H(r) is of 
course irrelevant, we retained a quadratic expression for convenience. In the case of the YOCP the bias function H(r) 
at (r, a* + 5a*) was chosen equal to the screening function H(r) obtained for the point (T, a*). 

To illustrate the method, we display in fig. (1) a plot of the bias function H(r) and of the unnormalized g w (r) 
at (r = 40, a* = 1). As a result of the structure of the bias function H(r) at short distance g w (r) exhibits a quite 
pronounced peak in this region. We show in fig. (2) the logarithm of the ratio K(r) of the unnormalized g u {r) by the 
properly normalized go(r) at the same point (r = 40, a* = 1). As can be seen in the figure the normalization constant 
of g u (r) can be obtained with a good accuracy by averaging K(r) for 1 < r/a < 2. The resulting function H(r) is 
displayed in fig. (3) for r varying in the range < r/a < 2. Note that in the region r m j„ < r/a < 2 where H(r) can be 
computed by a standard MC simulation it coincides almost perfectly with the H(r) computed by the biased scheme. 
In all cases the screening function can be accurately fitted by the simple polynomial ao + a2(r/a) 2 + a4(r/a) 4 + ae(r/a) e 
which traverses all the error bars. 

In the case of the OCP, the theoretical value a-i = — 1/4 8 " 10 is roughly recovered although its numerical value seems 
to depend strongly upon the number of particles. For instance at V — 40 one finds ai = —0.234, —0.239, —0.249, and 
—0.25 for samples of respectively N = 500, 1000, 2000, and 3000 particles. A more systematic study of these finite 
volume effects will be presented elsewhere. In the case of the OCP our results for H(r) differ by a small but significant 
amount of those of Ogata (see fig. (4)), but we have no firm conclusion concerning these discrepancies. A bunch of 
H(r) for various T is displayed in fig. (5). 

In the case of the YOCP the amplitude of the function H(r) at a given T decreases steadily as a increases as can 
be shown in fig. (6) but the shape of the curves at short distances remains unchanged. A 6 th order even polynomial 
perfectly fits the numerical data in the range < r/a < 2. A test of some theoretical predictions by Rosenfeld and 
Chabrier 11 concerning H(r) at small a is planned for future work. 
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FIG. 1. The bias function H(r) (top curve, solid line) used at (F = 40, a* = 1) to determine the function g w (r) (bottom 
curve). The error bars on g w (r) correspond to two standard deviations. 
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FIG. 2. Logarithm of the ratio K(r) of the normalized go(r) (standard MC calculation) and the unnormalized g u (r) (biased 
MC calculation) in the range r min < r < 2a at (r = 40, a* = 1). 
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FIG. 3. The screening function H(r) at (r = 40, a* = 1) for a sample of N = 500 particles. The standard MC calculation 
allows only the determination of H(r) in the range r min < r < 2a (curve on the right). The biased MC gives the curve in the 
range r < 2a. In the overlapping region the agreement between these two estimates is satisfactory. The solid curve is a 6 th 
order even polynomial of the MC data 
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FIG. 4. The screening function H(r) for the OCP at T = 40. Top curve: Ogata result. Bottom curves: MC data for 
sample of N = 500 ions (standard and biased simulations). The solid curve is a 6 th order even polynomial of the MC data. 
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